home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1995 April / macformat-023.iso / Shareware City / Developers / RLaB / toolbox / czt.r < prev    next >
Encoding:
Text File  |  1994-09-23  |  1.8 KB  |  91 lines  |  [TEXT/RLAB]
  1. //-------------------------------------------------------------------//
  2. // czt.r
  3.  
  4. // Syntax:    czt ( X , M , W , A )
  5.  
  6. // Description:
  7.  
  8. //    czt returns the M element Chirp z-transform, where:
  9. //    M: The number of points in the returned spectrum.
  10. //    W: The "spacing" between points along the z-pane spiral.
  11. //       A complex number along the unit circle.
  12. //    A: the complex starting point of the evaluated transform.
  13. //       A complex number along the unit circle.
  14.  
  15. //    At present X must be a row or column matrix. The column-wise
  16. //    matrix czt operation will be added later.
  17.  
  18. //      Ian Searle, 10/10/93
  19.  
  20. //    References:
  21. //    Rabiner, L.R. and B. Gold, Theory and Application of
  22. //     Digital Signal Processing, Prentice-Hall, Englewood
  23. //     Cliffs, New Jersey, pp. 393-399, 1975.
  24. //-------------------------------------------------------------------//
  25.  
  26. czt = function ( _x, M, W, A )
  27. {
  28.   global (pi)
  29.  
  30.   x = _x;
  31.   if (min (size (x)) != 1) { error ("czt: only vector X allowed"); }
  32.   nr = x.nr; nc = x.nc;
  33.  
  34.   if (nc == 1) 
  35.   { 
  36.     x = x'; 
  37.     nr = x.nr; nc = x.nc;
  38.     column = 1;
  39.   else
  40.     column = 0;
  41.   }
  42.   N = nc;
  43.  
  44.   if (!exist (M)) { M = length (x); }
  45.   if (!exist (W)) { W = exp (-1j .* 2 .* pi ./ M); }
  46.   if (!exist (A)) { A = 1; }
  47.  
  48.   //
  49.   // Compute a power of 2 value for fft()
  50.   //
  51.  
  52.   L = 1;
  53.   while (L < (N + M - 1)) { L = 2 .* L; }
  54.  
  55.   //
  56.   // Form the L-point sequence y(n)
  57.   //
  58.  
  59.   n = 0:(N-1);
  60.   y = [ A.^(-n) .* W.^((n.^2)./2) .* x , zeros (size (N+1:L)) ];
  61.  
  62.   //
  63.   // Form the L-point sequence v(n)
  64.   //
  65.  
  66.   n = 0:(M-1);
  67.   v = W.^(-(n.^2)./2);
  68.   v = [v, zeros (size (M+1:(L - N + 1)))];
  69.   n = (L-N+1):(L-1);
  70.   v[(L-N+2):L] = W.^(-((L-n).^2)./2);
  71.  
  72.   //
  73.   // Convolution
  74.   //
  75.  
  76.   Y = fft (y, L);
  77.   V = fft (v, L);
  78.   G = V.*Y;
  79.   g = ifft (G, L);
  80.  
  81.   k = 0:(M-1);
  82.   X = g[1:M] .* W.^((k.^2)./2);
  83.  
  84.   if (column)
  85.   {
  86.     return X';
  87.   else
  88.     return X;
  89.   }
  90. };
  91.